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Abstract 

We revisit here the mathematical model for ATP production in mitochondria 
introduced recently by Bertram, Pedersen, Luciani, and Sherman (BPLS) 
as a simplification of the more complete but intricate Magnus and Keizer's 
model. We correct some inaccuracies in the BPLS original approximations 
and then analyze some of the dynamical properties of the model. We infer 
from exhaustive numerical explorations that the enhanced BPLS equations 
have a unique attractor fixed point for physiologically acceptable ranges of 
mitochondrial variables and respiration inputs. We determine, in the station- 
ary regime, the dependence of the mitochondrial variables on the respiration 
inputs, namely the cytosolic concentration of calcium Cac and the substrate 
fructose 1,6-bisphosphate FBP. The same effect of calcium saturation re- 
ported for the original BPLS model is observed here. We find out, however, 
an interesting non- stationary effect: the inertia of the model tends to increase 
considerably for high concentrations of calcium. 

Keywords: Mitochondria, Calcium, ATP, Mathematical model 



1. Introduction 

The exchange of energy in cells is mostly mediated by ATP (adenosine 
triphosphate) molecules. Such molecules are produced in several processes in 
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an eukaryotic cell, but the principal source of ATP is typically the oxidative 
phosphorylation process which takes place in mitochondria. The mitochon- 
drion is an organelle with two membranes, having, therefore, two distinct bulk 
regions: the intermembrane space and the mitochondrial matrix. In the inner 
membrane, there are plenty of protein transporters and ionic channels, some 
of which execute an active transport leading to a gradient of some ions and 
molecules H]- The metabolic cascade that leads to the production of ATP 
in the mitochondrion starts in the cytoplasm. At first, glucose is transported 
from the extracellular medium into the cytoplasm by GLUT transporters. It 
is then converted in glucose-6-phosphate (G6P) by the enzyme hexokinase. 
G6P is then converted in pyruvate in a process called glycolysis, in which 
there is a net production of two ATP molecules. The pyruvate produced is 
transported into the mitochondrion (to the mitochondrial matrix) and is me- 
tabolized in a series of oxidation-reduction reactions in the citric acid cycle 
leading to the production of the nicotinamide adenine dinucleotide NAD and 
flavin adenine dinucleotide FAD. These electron donor molecules are oxidized 
in the complexes I to IV present in the inner mitochondrial membrane. These 
reactions lead to the activation of a proton pump, creating a pH gradient be- 
tween the inter membrane space and the matrix. The protons pumped into 
the intermembrane space return to the matrix through a transporter that 
uses their energy to catalyze the conversion of ADP (adenosine diphosphate) 
into ATP. The ATP produced in the mitochondria is then transported to the 
cytoplasm by the ATP / ADP exchanger [l], |2| . 

The kinetic aspects of the processes involved in the ATP production in 
mitochondria are rather intricate. This issue was addressed by Magnus and 
Keizer (MK), who introduced in the series of papers [s], 0, [sj a theoretical 
kinetic model for ATP production in mitochondria based on the known bio- 
physical properties of the enzymes and transporters involved in the process. 
In fact, the MK model was built by considering electrical activity and cy- 
tosolic calcium handling in insulin-secreting pancreatic /3-cells. The model 
consists basically in a set of equations describing the dynamics of the citric 
acid cycle, the proton pump, and the inner mitochondrial membrane trans- 
porters of ATP and calcium. The MK model is effectively based on first 
biophysical principles and provides a very detailed and accurate description 
of the processes considered to be important for mitochondrial oxidative phos- 
phorylation. However, it is also a rather complex model with cumbersome 
equations, preventing a systematic mathematical study of its dynamical and 
physiological properties. 
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A simplification of tlie MK model aiming to retain its main dynami- 
cal properties was introduced recently by Bertram, Pedersen, Luciani, and 
Sherman (BPLS) in [oj. The BPLS model incorporates some refinements in- 
troduced by Cortassa et al. in 0] for the description of the ATP production 
in cardiac cells. In fact, BPLS model can be considered as an approximation 
of the Cortassa et al.^s model instead of the original MK one. As we will see, 
this was probably the origin of some inaccuracies in the BPLS equations. As 
in the original MK model, the mitochondrial ATP production in the BPLS 
model is governed by four dynamical variables, namely the potential drop in 
the inner membrane AV^ and the mitochondrial concentrations of: reduced 
nicotinamide adenine dinucleotide NADH, adenosine diphosphate ADP, and 
calcium Ca^. The mitochondrial concentrations of pyridine and adenine nu- 
cleotides are assumed to be conserved 

NAD^ + NADH^ = NADtot, (1) 

ADP^ + ATP^ = Atot, (2) 

where NADtot and Atot stand for the total mitochondrial concentration of 
the respective nucleotides. The balance of the pertinent fluxes and reactions 
yields to the following equations 

^NADH = JpDH-Jo, (3) 



dt 

d 



where 



^^ADP = Jant-Jfifo, (4) 

— = fmi-Juni — JNaCa) , (5) 

= C*m^ i'^H — Ja'NT — <^NaCa " 2Juni) , (6) 



■Jh — ■JH,res — ■Jh,ATP — J H, leak- (7) 



The derivation and meaning of the fluxes presented in the right-handed sides 
of Eq. ©-([n]) are rather involved. The main details and the pertinent ref- 
erences can be found in the BPLS paper [6|. We have checked carefully the 
derivation of each of these fluxes and we have found out some inaccuracies 
in the BPLS expressions for adenine nucleotide translocator rate Jant and 
for calcium uniporter rate Juni- As we will see, these problems probably have 
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originated in the transcription of the original MK equations to the Cortassa 
et a/.'s model. 

In the present paper, we propose some enhanced approximations in the 
BPLS framework for the fluxes Jant and Juni and analyze some of the dynam- 
ical properties of Eqs. (IS])-©. We show, in particular, that for physiologically 
acceptable ranges of mitochondrial respiration inputs, namely the cytosolic 
concentration of calcium Cac and the substrate fructose 1,6-bisphosphate 
FBP, the BPLS equations have a unique physiologically acceptable attractor 
fixed point. Exhaustive numerical explorations indicate that the BPLS model 
is indeed globally stable, reinforcing its relevance to physiological quantita- 
tive studies, despite its simplicity when compared to the MK original model. 
We determine, in the stationary regime, the dependence on constant respira- 
tion inputs Cac and FBP of the of the four mitochondrial variables considered 
in the model. As in the original BPLS model, we observe here qualitatively 
distinct behavior for low and high Cac concentration excitations. We detect, 
moreover, a non-stationary effect: the inertia of the model tends to increase 
considerably for high concentrations of calcium. 



2. The Enhanced BPLS Model 

We will focus here in the problems for the BPLS expressions for the 
adenine nucleotide translocator rate Jant and for calcium uniporter rate 
Juni, since all the other quantities appearing in ([3])-® were checked to be 
correct and accurate for the physiological ranges of variables and parameters. 
The MK expression for the former is (see Eq. (16) of [sj) 

1 _ ATPc ADP^ 

o'ANT — V^max,ANT7 . x / n"- [o) 



The precise meaning of all the quantities presented in this formula can be 
found in js], 0, Isj or in the BPLS paper jsl as well. On the other hand, the 
expression for Jant presented in the Eq. (35) of the Cortassa et al. paper 
0] lacks the exponential in the numerator as in (jH]). The BPLS expression 
for Jant, obtained probably from Cortassa et al.'s Eq. (35), is 

•'ant-P^I^^U"-^^, (9) 

ADIV + ^20 / 
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where pig and p2o are numerical parameters. Notice that Eq. (35) of the 
BPLS paper has also a mistake in the denominator. The (reasonable) hy- 
pothesis used to simplify (|8]) in the BPLS model is the assumption that, due 
to the ion transporters action, the rates of ATP to ADP in the mitochondrial 
matrix and in the cytoplasm are approximately the same, 

ATPc ^ ATPjn (■]r\\ 
ADPe ADP^' ^ ' 

Note that this assumption implies from ([H]) that Jant ~ for ATPm — )■ Atot 
(and, hence, ADPm — )■ according to (E])), which is incompatible with the 
BPLS expression ([9]). The same occurs for the limit AV — )■ 0. With the 
assumption f lTOj) and taking into account the conservation of mitochondrial 
pyridine nucleotides ([2]), the adenine nucleotide translocator rate reads 

T -V ATP^ 1 - e-^ 

<>'ANT — »/max,ANT— ^ TTf^ Tfav ■ [^^) 

This is our first proposed approximation, which captures all the essential 
properties of ([H]) and is simple enough to be mathematically manipulated. 
Notice that for the typical range of physiological parameters, neglecting the 
exponential in the numerator of f|TT]) implies an relative error inferior to 5%. 
We will not, however, adopt this further approximation in this work. Figure 
([1]) illustrate the discrepancies between the expressions (Q and ffTTj) for typ- 
ical physiological ranges of the parameters and variables. An inspection of 
the graphics 9(B) of [6|] reveals that they have probably compared their ap- 
proximated expression with their Eq. (35), which was transcribed incorrectly 
from Cortassa's paper J3]- 

With respect to the calcium uniporter rate Juni, the original MK expres- 
sion reads (see Eq. (19) in j3|) 

_ j^jAV-AVo) + 



^-j^iAV-AVo) / _ ca, 



1 + + 



i^trans/ ' (1 + Cac /Xact 



In the BPLS derivation of the approximation for J^^i, it is used Eq. (38) of 
Cortassa et al. fi], where there is a mistake in the denominator. The BPLS 
proposed expression for the calcium uniporter rate is 

J,,i = (p2lA\/-p22)Cae^ (13) 
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Figure 1: Comparison between the orignal BPLS expression for the adenine nucleotide 
translocator rate Jant & and our proposal (fTTj) based on the MK model. Even the 
concavities of the curves are different. The curves were calculated by assuming ATP = 
3mM and Atot = 15 mM. An inspection of the graphics 9(B) of [6| reveals that they 
have probably compared their approximated expression with their Eq. (35), which was 
transcribed incorrectly from Cortassa's paper Q. 



where p2i and P22 are numerical parameters. We found this equation to be 
inaccurate for the typical physiological range of parameters. We propose to 
keep in the approximated model the complete original MK equation (1121) . Its 
dependence on AV is already in a rather simple form, and the complications 
for Cac are harmless for the dynamical analysis, as we will show. Now, by 
introducing the following dimensionless variables 

X = ^^PH- (14) 
ATP 



m 
tot 
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FBP , , 

(19) 



FBP, 



where the constant values are presented in Table [T|, and taking into account 
the new proposed expressions f|TT]) and (|T2l) . the rates in the right-handed 
sides of ([3])-® will be given by 

JpDU = riv^^^ (a2 + , (20) 

ai + z \ 1 — X J 

Jo = r2^—il + a,e''n-\ (21) 
a3 + X 

Jant = r^y- , (22) 

1-y 





= r,[{as + y){l + a,e"^^'^''")]'\ 


(23) 


J H,ves 


= CtllJo, 


(24) 


Jh,atp 


= ai2^FlF05 


(25) 


J H, leak 


= r5(w-ai3), 


(26) 






(27) 



u 

= (28) 

where 

= (29) 

ai8 + (l + ai6u)"-(l + ai7u)4 ^ ^ 

The values of the numerical parameters are presented in Table [H With the 
new dimensionless variables, the Eqs. ([3])- ([6]) can be cast in the form 

X = ( Jpi^H - Jo) , (30) 

(JpiFo — Jant) , (31) 



A 



iNADtot = 10 X lU /iiVl 


Atot = 15 X lU yUiVl 


Cao : 


= 0.2 /iiVl 


AVo = 91 mV 


FBPo = 1 /iM 


- 


= 1.8 /iMmV 


Kiax,ANT = 5 /iM ms 


Knax,uni = 10 yUM HIS 


F 
RT 


= 0.037 mV 


/ = 0.5 


-f^trans = 19/iM 


TV' 

-f^act 


= 0.38 /iM 


T 1 "1 n 

L = 110 


fm = 0.01 


na = 


2.8 


ri = U.z /iiVl ms 


r2 = 0.6 /iMms-^ 


r3 = 


5 /iM ms" 


r4 = 23.3;uMms"^ 


rs = 0.182 /iMms-i 


re = 


0.001 /iM ms-i 


r7 = 0.11 yuMms"^ 








ai = 0.05 


(22 = 1 


0-3 = 


0.01 


tti = 4.23 X 10-16 


ag = 18.2 


ag = 


3.37 


ar = 1.68 


ag = 0.67 


ag = 


5.10 X 10^ 


aio = 10.7 


ail = 11-7 


ai2 = 


= 3.43 


ai3 = 0.16 


ai4 = 1.46 


ai5 = 


= 6.73 


ai6 = 0.52 


ai7 = 0.01 


ai8 = 


= 110 



Table 1: Numerical parameters and rates for the enhanced BPLS model, see equations 
dUll-lall). AU the values were obtained from [3, Si] and 0. 



^ ~ ^ ^^NaCa) , (32) 

^ ~ AT/ ("^^ ~ <^ANT — <^NaCa " 2Juni) , (33) 

where /m and Cm stand, respectively, for the fraction of free Ca ions and 
the mitochondrial capacitance, see Table 1. Equations (!30!) - (l33l) form a non- 
autonomous systems of four first order differential equations. The external 
excitations u{t) and v{t) are related, respectively, to the cytosolic concen- 
tration of calcium Ca^ and the substrate fructose 1,6-bisphosphate FBP, see 
Eqs. fll8l) and (1191) . We can now start the dynamical analysis of the model. 



3. Dynamics of the model 

Let us consider initially the fixed points {x^,,y^, z^^w^) of the system (150]) - 
(155]) assuming constant inputs {u^^v^). The physiologically meaningful range 
for the variables x and y is [0, 1] by construction, see ([I])-© and (|T^ - (|T5]) . 
For z and w, we assume only that they are non negative. The typical phys- 
iological range for the potential drop, however, is more restrictive, corre- 
sponding to ^ [90,225]mV, which is equivalent to w ~ [1,2.5]. For 
the inputs u and f, we consider the ranges [0,10] and [0,20], respectively. 
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which corresponds to Cac ~ [0, 2] /xM and FBP [0, 20] /iM. We perform 
an exhaustive numerical search [sj for fixed points of fl30l)- fl33l) by assuming 
u G [0, 10] and v G [0, 20] constants. For all tested values of u and v, only 
one physiological {x, y G [0, 1] both z,w > 0) fixed point was found, which is 
always stable. Moreover, the fixed point is globally stable for physiological 
ranges of variables, meaning that any solution of fl5U]) - fl55]) with reasonable 
initial conditions will tend asymptotically to the fixed point. The variable w 
and X have the quickest convergence to the fixed point, where y and z are the 
slowest ones. The values of {x^,,y^, z^,w^) as function of the constant inputs 
{u^,v^:) are depicted in Fig. ([2]), from where one can already observe some 
physiologically consistent dynamical properties which we describe briefiy be- 
low. 

First, the production of ATP and the concentration of NADH vanishes 
in the absence of cytosolic calcium Cac and/or the substrate fructose 1,6- 
bisphosphate FBP, i.e., x and ?/ — )■ for m or t> — )■ 0. Also, we see that for 
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reasonable values of u and v the value of the potential drop AV (w) is almost 
constant and close to 150 mV. This stability is probably the reason why the 
original BPLS model is robust, despite the inaccuracies for the expression of 
Jant and Jo we are correcting in this paper. We will return to this point in 
the last section. 

Another important feature of the BPLS model is the reversion of some mi- 
tochondrial variables dynamical behavior in the presence of lower and higher 
concentration of calcium, which can be simulated, as described in [sj , by set- 
ting ai = in the expression for Jpdh (120|) . This conclusion can be inferred 
from Fig. |2l but it is better illustrated in Fig. [3l which depicts the solutions 
of (15n]) - (|55]) for an oscillatory Cac input of the form 

u{t) = uo + ui sin(t/to), (34) 

with constant v(t) and initial conditions {x{0),y{0), z{0),w{0)) given by the 
values of the fixed point corresponding to = m(0) and = "^(0). As we will 
see, such a choice of initial condition is consistent with the adiabatic regime 
we observe for sufficiently slow inputs (large period to). For lower values of 
uq (low Cac concentrations), all the mitochondrial variables increases and 
decreases in synchrony with the variations of u. On the other hand, for 
higher values of Uq, the dynamical behavior of x, y and w is reversed, i.e., 
they tend to decrease/increase while u increases/decreases. This effect can 
be understood from the relation between u^, and 2;* depicted in Fig. 2.c. The 
value of tends to increase rapidly when u^, increases, and this behavior can 
be traced back to the condition JNaCa = Juni defining the fixed point = 0. 
However, for large values of z, the dependence of the expression for Jpdh 
( 12U]) on z saturates and becomes equivalent to setting oi = 0. Without the z 
suppression term in Jpdh, the dynamical behavior of the variables x, y, and 
w is reversed, as it was pointed out in the original BPLS analysis. The value 
of V (FBP) does not change qualitatively this dynamical behavior. 

The oscillatory excitations illustrated in Fig. [3] has period to = 3 min. For 
inputs varying over a time scale of minutes, the system evolves adiabatically 
in a good approximation, i.e., the instantaneous solution {x(t),y{t), z{t),w{t)) 
is well approximated by the fixed point (x*, y^,z^:,w^) corresponding to -u* = 
u{t) and = v{t). In other words, for slowly varying inputs, the solutions 
of the system are confined to the fixed-point surfaces depicted in Fig. HJ 
Of course, one expects a breakdown of this adiabatic behavior for rapidly 
varying inputs. Non-stationary effects must appear for inputs varying with 
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1,63 




Figure 3: Response of the equations ([30 |) - ([33l) to oscillatory inputs (|34|) . Notice that for low 
Cac concentrations, all the mitochondrial variables increases and decreases in synchrony 
with the variations of u. On the other hand, for high Cac concentrations, the behavior of 
X, y and w is reversed. All the curves were evaluated for FBP = 0.5 fiM. See the text for 
further details. 
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a characteristic time smaller than a certain critical value. In order to study 
non-stationary effects, we consider the response of the system for inputs of 
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the type 

u{t) = Uq + Ui tanh(t/to), (35) 

for different values of to. This situation is depicted in Fig. IHfor some values 
of Uq and Ui and for to = 2.5, 1, 0.5, and 0.02 seconds. It is clear that for lower 
values of u (Cac), approximately 10 seconds are enough to assure that the 
system reaches an adiabatic regime. As we have already noticed, the variables 
y (ATP) and z (Cam) are the slowest ones to attain their respective stationary 
regime. Increasing the values of u implies the increasing of such "relaxation" 
time. The second column in Fig. H] corresponds to a situation with u G [2, 4], 
for which almost 30 seconds are necessary to assure the attainment of the 
stationary regime. The inertia of the system, hence, increases considerably 
for higher concentrations of cytosolic calcium. 

4. Final Remarks 

We have revisited here the mathematical model for ATP production in 
mitochondria introduced recently by Bertram, Pedersen, Luciani, and Sher- 
man (BPLS) in as a simplification of the more complete but intricate 
Magnus and Keizer's model [i], 0, [ij. We checked carefully all the approx- 
imations introduced in the BPLS model and found some inaccuracies for 
the approximations used for the adenine nucleotide translocator rate Jant 
and for calcium uniporter rate Juni- We proposed enhanced approximations 
for such rates based on the original Magnus and Keizer's model and ana- 
lyzed some dynamical properties of the model. Our results for the stationary 
regime indicate that the BPLS model is indeed globally stable, reinforcing 
its relevance to physiological quantitative studies, despite its simplicity when 
compared to the Magnus and Keizer's model. We have considered also the 
non-stationary regime and detected an interesting effect: the inertia of the 
model tends to increase considerably for high concentrations of cytosolic cal- 
cium. 

It is interesting to notice that the dynamics of our enhanced model are 
qualitatively similar to the original BPLS one, despite the differences in the 
rates Jant and Juni for physiological ranges, as depicted, for instance, in Fig. 
[U This point can be understood from the fact that the value of w*, which 
does not depends tightly on the details of such rates, is almost constant 
and corresponding to AV^ = 150 mV for reasonable values of the inputs 
and M*. For a fixed value of AV^, the numerical parameters in can be 
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fitted to provide a good adiustment for tlie real ATP dependence of (fTT]) . 
An inspection of Fig. 9 of [g] reveals that the adjustment of their numerical 
parameters was checked for ^ 160 mV, which is close to the physiological 
fixed point. 
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